#ifndef _REENTRANT
#define _REENTRANT
#endif // for using threads

/* STANDARD LIBRARY */
#include<cstdio>
#include<stdlib.h>
#include<iostream>
#include<string>
#include<math.h>
#include<pthread.h>
#include<errno.h>

/* ADDITIONAL LIBRARY */
#include"jk.h"
#include"constants.h" // constants
#include"supToGSL.h" // support to GSL library
#include"thread.h" // thread
#include"jkClass.h" // supprot to jk program (additional functions)
#include"jkFunc.h" // set of functions to solve using GSL
#include"colorText.h" // functions which change color of text (for fun)
using namespace std;

int main()
{
	unsigned int nThreads=1; // number of threads
	pthread_t *id[nThreads];
	pthread_mutex_t *mut = mutexNew();
	ParThre **parThre = new ParThre*[nThreads];
	parThre[0] = new ParThre();

	// STANDARD GA
	// double par1[7] = {0, 0.97, 1., 0.25, 500, 8, 1024};
	// double x01[3] = {-0.4, 0., 0.28};
	// STANDARD SGA
	// double par2[7] = {0, 0.97, 1., 0.25, 500, 8, 1024};
	// double x02[5] = {3.26, 0, 0.29, -3.65, 0.};


	// XV.3 2D, standard, testy rozdzielczosci: wykres od L;
	parThre[0]->MaxNumberOfShake = 5; // how many we can shake system?
	parThre[0]->maxInterp = 5; // how many points will be use to interpolation
	parThre[0]->mut = mut;
	parThre[0]->d = 6; // which parameter will be changed
	parThre[0]->nRepetitions = 32; // and how many times
	parThre[0]->start = 128; // from which value
	parThre[0]->end = 256; // to which value
	// SGA
	parThre[0]->nEquations = 5; // how many equations
	parThre[0]->func = &jk2dSGA;
	parThre[0]->title = "# 2D: SGA - testy rozdzielczosci - B\n# h\tstatus\titer\tmu\t\tm\t\td\t\tl_n\t\tl_m\n";
	const unsigned int nPar = 7; // 7 or 9
	const unsigned int nVar = 5;
	parThre[0]->nPar=nPar;
	double par[nThreads][nPar];
	double x[nThreads][nVar];

	unsigned int i=0;
	// 0
	parThre[i]->id = i; // id of current thread
	par[i] = {0.1, 0.97, 1., 0.25, 500, 8, 512};
	x[i] = {3.45, 0.351, 0.274, -3.68, 0.473};
	parThre[i]->par = par[i]; // parameters: h, n, t1, t2, beta, U, L
	parThre[i]->x0 = x[i]; // start point: // mu, m, d
	parThre[i]->file = "plots/2d-test-L/3dSGA_h01_L128-256_goodE.dat";
	id[i] = threadCreate(Run, parThre[i]);

//	// 1
//	i++;
//	parThre[i] = new ParThre(*parThre[i-1]);
//	// GA
//	parThre[i]->nEquations = 3; // how many equations
//	parThre[i]->func = &jk3dGA;
//	parThre[i]->title = "# 2D: GA - testy rozdzielczosci - B\n# h\tstatus\titer\tmu\t\tm\t\td\n";
//
//	parThre[i]->id = i; // id of current thread
//	par[i] = {0.1, 0.97, 1., 0.25, 500, 8, 1024};
//	x[i] = {-0.390, 0.0647, 0.287};
//	parThre[i]->par = par[i]; // parameters: h, n, t1, t2, beta, U, L
//	parThre[i]->x0 = x[i]; // start point: // mu, m, d
//	parThre[i]->file = "plots/2d-test-L/3dGA_h01_L128-256_goodE.dat";
//	id[i] = threadCreate(Run, parThre[i]);

	// Running
	mutexClose(mut);
	cleanScrean();
	textRed("Program MROS - test L dla 2 wymiarow, standard\n");
	mutexOpenUnlock(mut);

	for(unsigned int i=0; i<nThreads; i++)
	{
		//if(*id[i] == NULL) { continue; }
		fprintf(stderr, "[[%d:", i);
		threadJoinIgn(*id[i]);
		fprintf(stderr, "%d]]", i);
	}

	// Cleaning
	gotoxy(0,nThreads+2);
	mutexDestroy(mut);
	for(unsigned int i=0; i<nThreads; i++) { delete parThre[i]; }
	delete[] parThre;
	return 0;
}
//:~ end of main

void* Run(void *wsk)
{
	Clock *clock = new Clock();
	ParThre *parThre = static_cast<ParThre*>(wsk);
	int id = parThre->id;
	double *x0 = parThre->x0;
	double *par = parThre->par;
	unsigned int nPar = parThre->nPar;
	unsigned int d = parThre->d;
	unsigned int nRepetitions = parThre->nRepetitions;
	double start = parThre->start;
	double end = parThre->end;
	const char *file = static_cast<const char*>(parThre->file);
	const char *title = static_cast<const char*>(parThre->title);
	const size_t nEquations = parThre->nEquations;
	const unsigned int MaxNumberOfShake = parThre->MaxNumberOfShake;
	unsigned maxInterp = parThre->maxInterp;
	int (*func)(const gsl_vector*, void*, gsl_vector*) = parThre->func;
	pthread_mutex_t *mut = parThre->mut;

	FILE *fp = fopen(file, "w");
	if(fp==NULL)
	{
		fprintf(stderr, "Error with opening %s to write!\n", file);
		exit(EXIT_FAILURE);
	}

	Parameters Par(nPar, par); // list of parameters // 7 or 9!!!
	Parameters X0(nEquations, x0); // start point
	MultirootSolverF S(func, nEquations); 	// which function will be solved

	fprintf(fp, "# h=%.3f, n=%.3f, t1=%.3f, t2=%.3f", par[0], par[1], par[2], par[3]);
	fprintf(fp, " beta=%.3f, U=%.3f, L=%d\n", par[4], par[5], (int)par[6]);
	fprintf(fp, "%s", title);

	double **y = new double*[nEquations];
	for(unsigned int i=0; i<nEquations; i++) { y[i] = new double[nRepetitions+1]; }
	Parameters *Result = NULL;
	unsigned int nShake = 0;
	bool isStuck = 0;
	double *x = new double[nRepetitions+1];

	for(unsigned int i=0; i<=nRepetitions; i++)
	{
		// a critical section required - start
		mutexClose(mut);
		cleanLine(id+1);
		gotoxy(0,id+1);
		float bar = static_cast<float>(i)/static_cast<float>(nRepetitions);
		printf("|");
		for(unsigned int j=0; j<50; j++)
		{
			if(bar*50. <= static_cast<float>(j)) printf(" ");
			else printf("=");
		}
		printf("| ");
		int iswitch = i%4;
		switch(iswitch)
		{
			case 0: printf("|"); break;
			case 1: printf("/"); break;
			case 2: printf("-"); break;
			case 3: printf("\\"); break;
		}
		if(isStuck)
		{
			printf(" %.0f%%", bar*100);
			printf("  aborted :-(\n");
			mutexOpenUnlock(mut);
			for(unsigned int i=0; i<nEquations; i++) { delete []y[i]; }
			delete []y;
			delete clock;
			fclose(fp);
			return NULL; }
		if(bar < 1)
		{
			printf(" %.0f%%", bar*100);
			float time = clock->estimateEnd(bar);
			const char *descriptorT = (const char*)"sec";
			if(time>60) { time /= 60; descriptorT = (const char*)"min"; }
			else if(time>3600) { time /= 3600; descriptorT = (const char*)"hours"; }
			printf(" end in %.0f %s\n", time, descriptorT);
		}
		else { printf("  complete :-)\n"); }
		mutexOpenUnlock(mut);
		// a critical section required - end

		double h = (end-start)/static_cast<double>(nRepetitions);
		if(nRepetitions==0) h=0;
		Par[d] = static_cast<double>(i)*h + start;
		x[i] = Par[d];
		if(DEBUG) X0.print();
		S.solveFor(&Par, &X0);
		Result = S.getResult();

		if(Result == NULL) // shake of system
		{
			fprintf(stderr, "%d: NULL error for %.6f (%d time)\n", id, Par[d], nShake);
			X0.shake(++nShake);

			if(nShake>MaxNumberOfShake)
			{
				fprintf(stderr, "error for: Par[%d] = %.6f\n", d, Par[d]);
				fprintf(stderr, "It is imposible to solve the set of equations, sorry\n");
				isStuck = true;
			}

			i--;
			continue;
		}

		for(unsigned int j=0; j<nEquations; j++) { y[j][i] = (*Result)[j]; }

//		for(unsigned int j=0; j<nEquations; j++) { X0[j] = (*Result)[j]; }
		X0.interpolation(x, y, maxInterp, i, x[i], h);

		fprintf(fp, "%.9f\t", Par[d]);
		S.printResult(fp);
		nShake = 0;
		Result->deleteP(); delete Result;
	}

	for(unsigned int i=0; i<nEquations; i++) { delete []y[i]; }
	delete []y;
	delete clock;
	fclose(fp);
	return NULL;
}
